/*_glsl.h - The GLSL datatypes in C++
 *
 * Copyright(c) 2007 Bert van der Weerd<hypercharge@gmail.com>
 *
 * This is Free Software under BSD licence. Share and enjoy.
 */

/*
 * --------
 * OVERVIEW
 * 
 * This code does not depend on or use OpenGL nor GLSL in any way. it
 * just uses the GLSL API.
 *
 * Because the vec<T,N> and mat<T,N> types are templated to both
 * element type T and dimension N, the following is also possible:
 *
 *   vec< complex<float>, 200 >  myVector;
 *
 * Because both classes are implemented using std::tr1::array<T,N>,
 * they are both fully compatible with stl.
 *
 * ---------
 * BEHAVIOUR
 *
 * In pseudo-BNF-code, the 'vec' and 'mat' types can do this:
 * 
 * 1) Create
 *	Matrices or Vectors
 * of type
 *	[bool|int|float|double|complex<T>]
 * of dimension
 *	[N].
 *
 * 2) Provide GLSL api as [global|member] functions, where 'GLOBAL API' is the
 * original
 * 
 *	A = transpose(B);
 *	
 * syntax, and 'MEMBER API' is the
 * 
 *	A = B;
 *	A.transpose();
 *
 * syntax. The latter syntax can be used to avoid temporary object creation,
 * which can be too costly for bigger size (vectors|matrices). You can
 * mix both syntax forms.
 *
 * 3) Provide [(Read|Write) (Binary|Text)] input/output of the two
 * types. (Only text-ostream is implemented)
 *	
 * 4) These classes are written with these main goals in mind:
 *
 *	1) Easy API that uses the well-known GLSL syntax.
 *	2) Speed and performance by using CLEAN CODE.
 *	3) performance is optimized for vec3/mat3 type of user code.
 *	3) for bigger N, optimization and scaling are difficult.
 *
 * ------------
 * CODE EXAMPLE
 * 
 *	vec3 x = init_vec(1.0f,2.0f,3.0f), y;
 *	mat3 A,B,C;
 *
 *	C = diag(x);
 *	A = transpose(C);
 *	B = A + C;
 *
 *	y = (A * B * C) * x;
 *
 *	return length(y);
 *
 * -----------------------------
 * ADDITIONAL NON-GLSL FUNCTIONS
 * 
 *   diag(vec<T,N> v) - create diagonal matrix from vector
 *   diag(T x, int N) - create diagonal (NxN)-matrix from scalair
 *   transpose(vec<T,N,M>) - I dont see (yet) why this one is missing
 *                           from the original glsl spec.
 *
 * ------
 * ISSUES
 * 
 * There are things in the GL Shading Language that cannot be
 * expressed in C++. These are the Swizzle Operations and Component
 * Name Aliasses. Just Remember: This is C++, not GLSL.
 *
 * The only functions NOT implemented are the texture functions, the
 * fragment shader specific functions, and the noise functions.
 *
 * In gcc, the not() builtin cannot be used. There is an operator!()
 * that performs the same function and is available everywhere. Also,
 * the memberwise constructors are moved into init_mat() and init_vec()
 * functions because constructors cannot be partially instantiated.
 *
 * MATRICES ARE IMPLEMENTED IN ROW-MAJOR ORDER, as a single
 * array, which is the standard in OpenGL.
 *
 * --------------------
 * IMPLEMENTATION NOTES
 *
 * The code #includes only <cmath> and only calls it explicitly from
 * the std namespace, like in 'std::sin(M_PI)'. The #include <iostream>
 * is just for simple mat/vec text output.
 *
 * The code consists of the following main classes:
 *
 *	template <class T, int N> class vec { ... };
 *      template <class T, int N, int M=N> class mat { ... };
 *
 * This allows us to write:
 *
 *      typedef vec<float,3> vec3;
 *      typedef mat<float,4> mat4;
 *
 * ------------------
 * OPTIMIZATION NOTES
 *
 * If you go for the last byte and cycle, you'll object to the glsl
 * syntax that it creates temporary copies, for example in the function
 * below:
 * 
 *      inline vec<T,N> cos(const vec<T,N>& v)
 *      { vec<T,N> out; for (int i=0;i<N;i++) out[i] = cos(v[i]); return out; }
 *
 * Now, for those who want to avoid this, the whole glsl api is also
 * accesable as memberfunctions into vec<T,N> and mat<T,N>, like this:
 * 
 *      inline void cos(const vec<T,N>& v)
 *      { for (int i=0;i<N;i++) (*this)[i] = cos(v[i]); }
 *
 * -------------
 * WEB RESOURCES
 *
 * * Matrix Layouts -- http://www.mindcontrol.org/~hplus/graphics/matrix-layout.html
 * * C++ Custom Stream I/O -- http://www.parashift.com/c++-faq-lite/input-output.html
 *
 * ----
 * TODO
 * 
 * 1) bring mat<> up to par with saru imp
 * 2) compare ColorRGB[A] with glsl api and bring up to par
 * 3) Read & up-to-par with http://www.morrowland.com/apron/tutorials/gl/gl_matrix.php
 *
 * -----
 * IDEAS
 * 
 * - stream reading of vec/mat
 * - combine stream read/write with ColMaj matrices so you can
 * - boolean vectors and logic for matrices and optimize using saru
 *   principles.
 * - implement Column Major matrices by partial templ. specialization.
 * - try out these classes a pure .NET code by vs2005 and check with
 *   saru for differences.
 * - std::vector integration might be a good idea? (low-priotiy).
 *
 */

#ifndef INCLUDED_GLSL_H
#define INCLUDED_GLSL_H


#include <iondrive/math.h>

#include <iostream> // so we can write to stdout streams

#ifndef INCLUDED_GLSL_NO_ARRAY
# include <tr1/array>
#endif


namespace iondrive {
namespace glsl {

	using namespace iondrive;



/************************************************************
 *
 * scalar builtin functions
 * 
 */





	//
	// Scalair functions defined by GLSL
	// 


	template <class T> inline T radians(const T& degrees)                                         { return T(iondrive::math::m_pi/180)*degrees; }
	template <class T> inline T degrees(const T& radians)                                         { return T(180/iondrive::math::m_pi)*radians; }
	template <class T> inline T exp2(const T& x)                                                  { return T(iondrive::math::exp(x * iondrive::math::m_ln2)); }
	template <class T> inline T log2(const T& x)                                                  { return T(iondrive::math::log(x) / iondrive::math::m_ln2); }
	template <class T> inline T inversesqrt(const T& x)                                           { return 1/iondrive::math::sqrt(x); }
	template <class T> inline T sign(const T& x)                                                  { return T((x>0) ? T(1) : ((x<0) ? T(-1):T(0))); }
	template <class T> inline T fract(const T& x)                                                 { return x - iondrive::math::floor(x); }
	template <class T> inline T mod(const T& x, const T& y)                                       { return T(iondrive::math::fmod(x, y)); }
	template <class T> inline T min(const T& x, const T& y)                                       { return y < x ? y : x; }
	template <class T> inline T max(const T& x, const T& y)                                       { return x < y ? y : x; }
	template <class T> inline T clamp(const T& x, const T& minVal = T(0), const T& maxVal = T(1)) { return glsl::min(glsl::max(x,minVal),maxVal); }
	template <class T> inline T mix(const T& x, const T& y, const T& a)                           { return x*(1-a) + y*a; }
	template <class T> inline T step(const T& edge, const T& x)                                   { return x<=edge ? T(0) : T(1); }
	template <class T> inline T smoothstep(const T& edge0,const T& edge1,const T& x)              { T t = clamp((x-edge0) / (edge1-edge0), T(0), T(1)); return t * t * (3 - 2*t); }
	template <class T> inline T length(const T& x)                                                { return iondrive::math::sqrt(x*x); }
	template <class T> inline T distance(const T& p0, const T& p1)                                { return length(p0-p1); }
	template <class T> inline T dot(const T& x, const T& y)                                       { return x*y; }
	template <class T> inline T normalize(const T&)                                               { return T(1); } // this is not the most useful function in the world
	template <class T> inline T faceforward(const T& N, const T& l, const T& Nref)                { return Nref*l < 0 ? N : -N; }
	template <class T> inline T reflect(const T& I, const T& N)                                   { return I - T(2)*N*I*N; }







/************************************************************
 *
 * Sequence functions
 * 
 */

template <class T, int N> struct vec;	// Forward declaration

// Functions are mostly the same, so a #define macro saves typing here.

#define SEQUENCE_FUNC01(__func_name)                                                    \
  template <class InputIterator,class OutputIterator>                                   \
  inline void __func_name(InputIterator begin, InputIterator end, OutputIterator out)   \
    {                                                                                   \
      while (begin!=end) *out++ = __func_name(*begin++);                                \
    }                                                                                   \
  template <class T, int N>                                                             \
  inline vec<T,N> __func_name(const vec<T,N>& v)                                        \
    {                                                                                   \
      vec<T,N> out;  /* this temp may or may not be optimized away */                   \
      __func_name(v.begin(),v.end(), out.begin());                                      \
      return out;                                                                       \
    }

#define SEQUENCE_FUNC02(__func_name,__func_fqn)                                         \
  template <class InputIterator,class OutputIterator>                                   \
  inline void __func_name(InputIterator begin, InputIterator end, OutputIterator out)   \
    {                                                                                   \
      while (begin!=end) *out++ = __func_fqn(*begin++);                                 \
    }                                                                                   \
  template <class T, int N>                                                             \
  inline vec<T,N> __func_name(const vec<T,N>& v)                                        \
    {                                                                                   \
      vec<T,N> out;                                                                     \
      __func_name(v.begin(),v.end(), out.begin());                                      \
      return out;                                                                       \
    }

#define SEQUENCE_FUNC03(__func_name,__func_fqn)                                         \
  template <class InputIterator,class OutputIterator>                                   \
  inline void __func_name(InputIterator a_begin, InputIterator a_end,                   \
                          InputIterator b_begin,                                        \
			  OutputIterator out)                                           \
    {                                                                                   \
      while (a_begin!=a_end) *out++ = __func_fqn(*a_begin++,*b_begin++);                \
    }

#define SEQUENCE_FUNC04(__func_name,__func_fqn)                                         \
  template <class InputIterator,class OutputIterator, class T>                          \
  inline void __func_name(InputIterator a_begin, InputIterator a_end,                   \
	                  const T& t,                                                   \
			  OutputIterator out)                                           \
    {                                                                                   \
      while (a_begin!=a_end) *out++ = __func_fqn(*a_begin++,t);                         \
    }                                                                                   \
  SEQUENCE_FUNC03(__func_name,__func_fqn)

SEQUENCE_FUNC01(radians)
SEQUENCE_FUNC01(degrees)
SEQUENCE_FUNC02(sin, iondrive::math::sin)
SEQUENCE_FUNC02(cos, iondrive::math::cos)
SEQUENCE_FUNC02(tan, iondrive::math::tan)
SEQUENCE_FUNC02(asin, iondrive::math::asin)
SEQUENCE_FUNC02(acos, iondrive::math::acos)
SEQUENCE_FUNC02(atan, iondrive::math::atan)
SEQUENCE_FUNC03(atan2, iondrive::math::atan2)
SEQUENCE_FUNC04(pow, iondrive::math::pow)
SEQUENCE_FUNC01(exp2)
SEQUENCE_FUNC01(log2)
SEQUENCE_FUNC02(exp, iondrive::math::exp)
SEQUENCE_FUNC02(log, iondrive::math::log)
SEQUENCE_FUNC02(sqrt, iondrive::math::sqrt)
SEQUENCE_FUNC01(inversesqrt)
SEQUENCE_FUNC02(abs, iondrive::math::fabs)
SEQUENCE_FUNC01(sign)
SEQUENCE_FUNC02(floor, iondrive::math::floor)
SEQUENCE_FUNC02(ceil, iondrive::math::ceil)
SEQUENCE_FUNC01(fract)
SEQUENCE_FUNC04(mod, glsl::mod)
SEQUENCE_FUNC04(min, glsl::min)
SEQUENCE_FUNC04(max, glsl::max)



/************************************************************
 *
 * vector class
 * 
 */



template <class T, int N>
#ifndef INCLUDED_GLSL_NO_ARRAY
class vec : public std::tr1::array<T,N>
#else
class vec
#endif
{
	
#ifndef INCLUDED_GLSL_NO_ARRAY

	// std::tr1::array<> handles the data storage of the array etc.
	
#else

	private:
		T data[N];
		
	public:
		// *all* access is done through this operator
		inline T& operator[](const int i)       { return data[i]; }
		inline const T&  operator[](const int i) const { return data[i]; }

#endif
		
	public:

		//
		// constructors
		//
		
		inline vec()                         { for (int i=0;i<N;i++)   (*this)[i] = T(); }
		inline vec(const T& t)               { for (int i=0;i<N;i++)   (*this)[i] = t; }
		inline vec(vec<T,N-1> a, const T& t) { for (int i=0;i<N-1;i++) (*this)[i] = a[i];  (*this)[N-1] = t; } // dimension up cast
		inline vec(vec<T,N+1> a)             { for (int i=0;i<N;i++)   (*this)[i] = a[i];  }  // dimension down cast
		
		
		//
		// experimental : non-temp-copy-creating interface
		//
		//
		// these inline functions below contain the glsl
		// functionality in member function form, which allows
		// you to use the vector class without generating the
		// temporary copies that inevitably arise when using
		// the glsl functional expression interface.
		//
		// ------------
		// CODE EXAMPLE
		//
		// vec3 x,y,z;
		// mat3 A,B;
		//
		// x = sin(A*y) + cos(B*z);
		//
		// ==[ BECOMES ]=======================
		// vec3 tmp1,tmp2,tmp3;
		//
		// tmp1.mul(A,y);                      // not implemented yet :(
		// tmp2.sin(tmp1);
		// tmp1.mul(B,z);
		// tmp3.cos(tmp1);
		// x.add(tmp2,tmp3);
		//
		//
		// 

		
		
		// init

		inline void init_vec(const T arg) { for (int i=0;i<N;i++) (*this)[i] = arg; }
		inline void init_vec(const T arg[]) { for (int i=0;i<N;i++) (*this)[i] = arg[i]; }
		inline void init_vec(const vec<T,N>& arg) { for (int i=0;i<N;i++) (*this)[i] = arg[i]; }

		
		// operators (result,v,a,b are vectors, t is scalair)

		
		// result.xxx(v)
		inline void add(const vec<T,N>& rhs) { for (int i=0;i<N;i++) (*this)[i] += rhs[i]; }
		inline void sub(const vec<T,N>& rhs) { for (int i=0;i<N;i++) (*this)[i] -= rhs[i]; }
		inline void mul(const vec<T,N>& rhs) { for (int i=0;i<N;i++) (*this)[i] *= rhs[i]; }
		inline void div(const vec<T,N>& rhs) { for (int i=0;i<N;i++) (*this)[i] /= rhs[i]; }
		// result.xxx(a,b);
		inline void add(const vec<T,N>& lhs,const vec<T,N>& rhs) { for (int i=0;i<N;i++) (*this)[i] = lhs[i] + rhs[i]; }
		inline void sub(const vec<T,N>& lhs,const vec<T,N>& rhs) { for (int i=0;i<N;i++) (*this)[i] = lhs[i] - rhs[i]; }
		inline void mul(const vec<T,N>& lhs,const vec<T,N>& rhs) { for (int i=0;i<N;i++) (*this)[i] = lhs[i] * rhs[i]; }
		inline void div(const vec<T,N>& lhs,const vec<T,N>& rhs) { for (int i=0;i<N;i++) (*this)[i] = lhs[i] / rhs[i]; }
		// result.xxx(t);
		inline void add(const T& rhs) { for (int i=0;i<N;i++) (*this)[i] += rhs; }
		inline void sub(const T& rhs) { for (int i=0;i<N;i++) (*this)[i] -= rhs; }
		inline void mul(const T& rhs) { for (int i=0;i<N;i++) (*this)[i] *= rhs; }
		inline void div(const T& rhs) { for (int i=0;i<N;i++) (*this)[i] /= rhs; }
		// result.xxx(v,t);
		inline void add(const vec<T,N>& lhs,const T& rhs) { for (int i=0;i<N;i++) (*this)[i] = lhs[i] + rhs; }
		inline void sub(const vec<T,N>& lhs,const T& rhs) { for (int i=0;i<N;i++) (*this)[i] = lhs[i] - rhs; }
		inline void mul(const vec<T,N>& lhs,const T& rhs) { for (int i=0;i<N;i++) (*this)[i] = lhs[i] * rhs; }
		inline void div(const vec<T,N>& lhs,const T& rhs) { for (int i=0;i<N;i++) (*this)[i] = lhs[i] / rhs; }

		inline void neg() { for (int i=0;i<N;i++) (*this)[i] = - ((*this)[i]); } // result.neg();
		inline void neg(const vec<T,N>& v) { for (int i=0;i<N;i++) (*this)[i] = - (v[i]); } // result.neg(vec_to_negate);


		// glsl functions
		

		inline void radians()                   { glsl::radians(this->begin(),this->end(), this->begin()); }
		inline void radians(const vec<T,N>& v)  { glsl::radians(v.begin(),v.end(), this->begin()); }
		
 		inline void degrees()                                                                   { for (int i=0;i<N;i++) (*this)[i] = iondrive::glsl::degrees((*this)[i]); }
		inline void degrees(const vec<T,N>& v)                                                  { for (int i=0;i<N;i++) (*this)[i] = iondrive::glsl::degrees(v[i]); }
		inline void sin()                                                                       { for (int i=0;i<N;i++) (*this)[i] = iondrive::math::sin((*this)[i]); }
		inline void sin(const vec<T,N>& v)                                                      { for (int i=0;i<N;i++) (*this)[i] = iondrive::math::sin(v[i]); }
		inline void cos()                                                                       { for (int i=0;i<N;i++) (*this)[i] = iondrive::math::cos((*this)[i]); }
		inline void cos(const vec<T,N>& v)                                                      { for (int i=0;i<N;i++) (*this)[i] = iondrive::math::cos(v[i]); }
		inline void tan()                                                                       { for (int i=0;i<N;i++) (*this)[i] = iondrive::math::tan((*this)[i]); }
		inline void tan(const vec<T,N>& v)                                                      { for (int i=0;i<N;i++) (*this)[i] = iondrive::math::tan(v[i]); }
		inline void asin()                                                                      { for (int i=0;i<N;i++) (*this)[i] = iondrive::math::asin((*this)[i]); }
		inline void asin(const vec<T,N>& v)                                                     { for (int i=0;i<N;i++) (*this)[i] = iondrive::math::asin(v[i]); }
		inline void acos()                                                                      { for (int i=0;i<N;i++) (*this)[i] = iondrive::math::acos((*this)[i]); }
		inline void acos(const vec<T,N>& v)                                                     { for (int i=0;i<N;i++) (*this)[i] = iondrive::math::acos(v[i]); }
		inline void atan()                                                                      { for (int i=0;i<N;i++) (*this)[i] = iondrive::math::atan((*this)[i]); }
		inline void atan(const vec<T,N>& v)                                                     { for (int i=0;i<N;i++) (*this)[i] = iondrive::math::atan(v[i]); }
		inline void atan2(const vec<T,N>& rhs)                                                  { for (int i=0;i<N;i++) (*this)[i] = iondrive::math::atan2((*this)[i],rhs[i]); }
		inline void atan(const vec<T,N>& v,const vec<T,N>& u)                                   { for (int i=0;i<N;i++) (*this)[i] = iondrive::math::atan2(v[i],u[i]); }
		inline void pow(const vec<T,N>& rhs)                                                    { for (int i=0;i<N;i++) (*this)[i] = iondrive::math::pow((*this)[i],rhs[i]); }
		inline void pow(const vec<T,N>& lhs,const vec<T,N> rhs)                                 { for (int i=0;i<N;i++) (*this)[i] = iondrive::math::pow(lhs[i],rhs[i]); }
		inline void exp2()                                                                      { for (int i=0;i<N;i++) (*this)[i] = iondrive::glsl::exp2((*this)[i]); }
		inline void exp2(const vec<T,N>& v)                                                     { for (int i=0;i<N;i++) (*this)[i] = iondrive::glsl::exp2(v[i]); }
		inline void log2()                                                                      { for (int i=0;i<N;i++) (*this)[i] = iondrive::glsl::log2((*this)[i]); }
		inline void log2(const vec<T,N>& v)                                                     { for (int i=0;i<N;i++) (*this)[i] = iondrive::glsl::log2(v[i]); }
		inline void exp()                                                                       { for (int i=0;i<N;i++) (*this)[i] = iondrive::math::exp((*this)[i]); }
		inline void exp(const vec<T,N>& v)                                                      { for (int i=0;i<N;i++) (*this)[i] = iondrive::math::exp(v[i]); }
		inline void log()                                                                       { for (int i=0;i<N;i++) (*this)[i] = iondrive::math::log((*this)[i]); }
		inline void log(const vec<T,N>& v)                                                      { for (int i=0;i<N;i++) (*this)[i] = iondrive::math::log(v[i]); }
		inline void sqrt()                                                                      { for (int i=0;i<N;i++) (*this)[i] = iondrive::math::sqrt((*this)[i]); }
		inline void sqrt(const vec<T,N>& v)                                                     { for (int i=0;i<N;i++) (*this)[i] = iondrive::math::sqrt(v[i]); }
		inline void inversesqrt()                                                               { for (int i=0;i<N;i++) (*this)[i] = iondrive::glsl::inversesqrt((*this)[i]); }
		inline void inversesqrt(const vec<T,N>& v)                                              { for (int i=0;i<N;i++) (*this)[i] = iondrive::glsl::inversesqrt(v[i]); }
		inline void abs()                                                                       { for (int i=0;i<N;i++) (*this)[i] = iondrive::math::abs((*this)[i]); }
		inline void abs(const vec<T,N>& v)                                                      { for (int i=0;i<N;i++) (*this)[i] = iondrive::math::abs(v[i]); }
		inline void sign()                                                                      { for (int i=0;i<N;i++) (*this)[i] = iondrive::glsl::sign((*this)[i]); }
		inline void sign(const vec<T,N>& v)                                                     { for (int i=0;i<N;i++) (*this)[i] = iondrive::glsl::sign(v[i]); }
		inline void floor()                                                                     { for (int i=0;i<N;i++) (*this)[i] = iondrive::math::floor((*this)[i]); }
		inline void floor(const vec<T,N>& v)                                                    { for (int i=0;i<N;i++) (*this)[i] = iondrive::math::floor(v[i]); }
		inline void ceil()                                                                      { for (int i=0;i<N;i++) (*this)[i] = iondrive::math::ceil((*this)[i]); }
		inline void ceil(const vec<T,N>& v)                                                     { for (int i=0;i<N;i++) (*this)[i] = iondrive::math::ceil(v[i]); }
		inline void fract()                                                                     { for (int i=0;i<N;i++) (*this)[i] = iondrive::glsl::fract((*this)[i]); }
		inline void fract(const vec<T,N>& v)                                                    { for (int i=0;i<N;i++) (*this)[i] = iondrive::glsl::fract(v[i]); }
		inline void mod(const T& rhs)                                                           { for (int i=0;i<N;i++) (*this)[i] = iondrive::glsl::mod((*this)[i],rhs); }
		inline void mod(const vec<T,N>& lhs,const vec<T,N>& rhs)                                { for (int i=0;i<N;i++) (*this)[i] = iondrive::glsl::mod(lhs[i],rhs[i]); }
		inline void min(const T& rhs)                                                           { for (int i=0;i<N;i++) (*this)[i] = iondrive::glsl::min((*this)[i],rhs); }
		inline void min(const vec<T,N>& lhs,const vec<T,N>& rhs)                                { for (int i=0;i<N;i++) (*this)[i] = iondrive::glsl::min(lhs[i],rhs[i]); }
		inline void max(const T& rhs)                                                           { for (int i=0;i<N;i++) (*this)[i] = iondrive::glsl::max((*this)[i],rhs); }
		inline void max(const vec<T,N>& lhs,const vec<T,N>& rhs)                                { for (int i=0;i<N;i++) (*this)[i] = iondrive::glsl::max(lhs[i],rhs[i]); }
		inline void clamp(const T& minVal = T(0), const T& maxVal = T(1))                       { for (int i=0;i<N;i++) (*this)[i] = iondrive::glsl::clamp((*this)[i],minVal,maxVal); }
		inline void clamp(const vec<T,N>& v, const T& minVal, const T& maxVal)                  { for (int i=0;i<N;i++) (*this)[i] = iondrive::glsl::clamp(v[i],minVal,maxVal); }
		inline void clamp(const vec<T,N>& minVal, const vec<T,N>& maxVal)                       { for (int i=0;i<N;i++) (*this)[i] = iondrive::glsl::clamp((*this)[i], minVal[i], maxVal[i]); }
		inline void clamp(const vec<T,N>& v,const vec<T,N>& minVal, const vec<T,N>& maxVal)     { for (int i=0;i<N;i++) (*this)[i] = iondrive::glsl::clamp(v[i],minVal[i],maxVal[i]); }
		inline void mix(const vec<T,N>& a, const vec<T,N>& b, const T& t)                       { for (int i=0;i<N;i++) (*this)[i] = iondrive::glsl::mix(a[i],b[i],t); }
		inline void mix(const vec<T,N>& a, const vec<T,N>& b, const vec<T,N>& c)                { for (int i=0;i<N;i++) (*this)[i] = iondrive::glsl::mix(a[i],b[i],c[i]); }
		inline void step(const vec<T,N>& edge, const vec<T,N>& x)                               { for (int i=0;i<N;i++) (*this)[i] = iondrive::glsl::step(edge[i],x[i]); }
		inline void smoothstep(const vec<T,N>& edge0, const vec<T,N>& edge1, const vec<T,N>& x) { for (int i=0;i<N;i++) (*this)[i] = iondrive::glsl::smoothstep(edge0[i],edge1[i],x[i]); }

		inline T length() const
		{
			return T(iondrive::math::sqrt(dot(*this)));
		}
		inline void normalize()
		{
			T len = length();
			if (len == 0) return;
			div(len);
		}
		inline void normalize(const vec<T,N>& v)
		{
			T len = v.length();
			if (len == 0) return;
			div(v,len);
		}

		inline void faceforward(const vec<T,N>& n, const vec<T,N>& i, const vec<T,N>& nref)
		{
			// return dot(nref,i) < 0 ? n : -n;

			T tmp1 = nref.dot(i);
			if (tmp1 < 0)
				for(int tmp2=0;tmp2<N;tmp2++) (*this)[tmp2] = n[tmp2];
			else
				for(int tmp3=0;tmp3<N;tmp3++) (*this)[tmp3] = -n[tmp3];
		}
		inline void reflect(const vec<T,N>& i, const vec<T,N>& n)
		{
			// return i - 2*dot(n,i)*n;
			vec<T,N> tmp2;
			tmp2.mul(n,2*n.dot(i));
			this->sub(i,tmp2);
		}

		inline T distance(const vec<T,N>& a, const vec<T,N>& b)
		{
			// return length(a-b);
			vec<T,N> tmp1;
			tmp1.sub(a,b);
			return T(tmp1.length());
		}
		inline T dot(const vec<T,N>& rhs) const
		{
			// T t; for (int i = 0; i < N; i++) t += a[i]*b[i]; return t;

			T t;
			for (int i=0;i<N;i++) t += (*this)[i]*rhs[i];
			return t;
		}
		inline void cross(const vec<T,3>& a,const vec<T,3>& b)
		{
			(*this)[0] = a[1]*b[2] - a[2]*b[1];
			(*this)[1] = a[2]*b[0] - a[0]*b[2];
			(*this)[2] = a[0]*b[1] - a[1]*b[0];
		}
		inline void cross(const vec<T,3>& b) // above cross() could be more efficient
		{
			vec<T,3> a;

			a[0] = (*this)[0];
			a[1] = (*this)[1];
			a[2] = (*this)[2];

			(*this)[0] = a[1]*b[2] - a[2]*b[1];
			(*this)[1] = a[2]*b[0] - a[0]*b[2];
			(*this)[2] = a[0]*b[1] - a[1]*b[0];
		}

		// Additional saru functions
		inline T squared_length() { T tmp; for (int i=0;i<N;i++) tmp += (*this)[i]*(*this)[i]; return tmp; }
		inline vec<T,3> unit_cross(const vec<T,3>& b)
		{
			vec<T,3> result;

			result[0] = (*this)[1]*b[2] - (*this)[2]*b[1];
			result[1] = (*this)[2]*b[0] - (*this)[0]*b[2];
			result[2] = (*this)[0]*b[1] - (*this)[1]*b[0];

			result.normalize();

			return result;
		}
};


// C++ LIMIT: We cannot do partial template specialization in the
// constructor, so these init_vec[234]() functions are provided as a fix.

template <class T> inline vec<T,2> init_vec(const T& x, const T& y)                         { vec<T,2> out; out[0] = x; out[1] = y; return out; }
template <class T> inline vec<T,3> init_vec(const T& x, const T& y, const T& z)             { vec<T,3> out; out[0] = x; out[1] = y; out[2] = z; return out; }
template <class T> inline vec<T,4> init_vec(const T& x, const T& y, const T& z, const T& w) { vec<T,4> out; out[0] = x; out[1] = y; out[2] = z; out[3] = w; return out; }

// REQUIRED: Vector operators

template <class T, int N> inline vec<T,N> operator+(const vec<T,N>& a, const vec<T,N>& b) { vec<T,N> out; for (int i=0;i<N;i++) out[i] = a[i]+b[i]; return out; }
template <class T, int N> inline vec<T,N> operator-(const vec<T,N>& a, const vec<T,N>& b) { vec<T,N> out; for (int i=0;i<N;i++) out[i] = a[i]-b[i]; return out; }
template <class T, int N> inline vec<T,N> operator*(const vec<T,N>& a, const vec<T,N>& b) { vec<T,N> out; for (int i=0;i<N;i++) out[i] = a[i]*b[i]; return out; }
template <class T, int N> inline vec<T,N> operator/(const vec<T,N>& a, const vec<T,N>& b) { vec<T,N> out; for (int i=0;i<N;i++) out[i] = a[i]/b[i]; return out; }

template <class T, int N> inline vec<T,N> operator*(const T& t, const vec<T,N>& v)        { vec<T,N> out; for (int i=0;i<N;i++) out[i] = t*v[i]; return out; }
template <class T, int N> inline vec<T,N> operator*(const vec<T,N>& v, const T& t)        { vec<T,N> out; for (int i=0;i<N;i++) out[i] = v[i]*t; return out; }
template <class T, int N> inline vec<T,N> operator+(const T& t, const vec<T,N>& v)        { vec<T,N> out; for (int i=0;i<N;i++) out[i] = t+v[i]; return out; }
template <class T, int N> inline vec<T,N> operator+(const vec<T,N>& v, const T& t)        { vec<T,N> out; for (int i=0;i<N;i++) out[i] = v[i]+t; return out; }
template <class T, int N> inline vec<T,N> operator-(const T& t, const vec<T,N>& v)        { vec<T,N> out; for (int i=0;i<N;i++) out[i] = t-v[i]; return out; }
template <class T, int N> inline vec<T,N> operator-(const vec<T,N>& v, const T& t)        { vec<T,N> out; for (int i=0;i<N;i++) out[i] = v[i]-t; return out; }
template <class T, int N> inline vec<T,N> operator/(const vec<T,N>& v, const T& t)        { vec<T,N> out; T tmp=T(1)/t; for (int i=0;i<N;i++) out[i] = v[i]*tmp; return out; } // speedup from saru

template <class T, int N> inline vec<T,N>& operator+=(vec<T,N>& a, const vec<T,N>& b)     { a = a + b; return a; } // this looks improvable
template <class T, int N> inline vec<T,N>& operator-=(vec<T,N>& a, const vec<T,N>& b)     { a = a - b; return a; }
template <class T, int N> inline vec<T,N>& operator*=(vec<T,N>& a, const vec<T,N>& b)     { a = a * b; return a; }
template <class T, int N> inline vec<T,N>& operator/=(vec<T,N>& a, const vec<T,N>& b)     { a = a / b; return a; }
template <class T, int N> inline vec<T,N>& operator+=(vec<T,N>& a, const T& t)            { a = a + t; return a; }
template <class T, int N> inline vec<T,N>& operator-=(vec<T,N>& a, const T& t)            { a = a - t; return a; }
template <class T, int N> inline vec<T,N>& operator*=(vec<T,N>& a, const T& t)            { a = a * t; return a; }
template <class T, int N> inline vec<T,N>& operator/=(vec<T,N>& a, const T& t)            { a = a / t; return a; }
// unary minus operator
template <class T, int N> inline vec<T,N> operator-(const vec<T,N>& v)                    { vec<T,N> out; for (int i = 0; i < N; i++) out[i] = -v[i]; return out; }

// REQUIRED: Builtin GLSL functions



//template <class T, int N> inline vec<T,N> radians(const vec<T,N>& v) { vec<T,N> out; radians(v.begin(),v.end(), out.begin()); return out; }
//template <class T, int N> inline vec<T,N> degrees(const vec<T,N>& v)                                                  { vec<T,N> out; for (int i = 0; i < N; i++) out[i] = degrees(v[i]); return out; }
//template <class T, int N> inline vec<T,N> sin(const vec<T,N>& v)                                                      { vec<T,N> out; for (int i = 0; i < N; i++) out[i] = iondrive::math::sin(v[i]); return out; }
//template <class T, int N> inline vec<T,N> cos(const vec<T,N>& v)                                                      { vec<T,N> out; for (int i = 0; i < N; i++) out[i] = iondrive::math::cos(v[i]); return out; }
//template <class T, int N> inline vec<T,N> tan(const vec<T,N>& v)                                                      { vec<T,N> out; for (int i = 0; i < N; i++) out[i] = iondrive::math::tan(v[i]); return out; }
//template <class T, int N> inline vec<T,N> asin(const vec<T,N>& v)                                                     { vec<T,N> out; for (int i = 0; i < N; i++) out[i] = iondrive::math::asin(v[i]); return out; }
//template <class T, int N> inline vec<T,N> acos(const vec<T,N>& v)                                                     { vec<T,N> out; for (int i = 0; i < N; i++) out[i] = iondrive::math::acos(v[i]); return out; }
//template <class T, int N> inline vec<T,N> atan(const vec<T,N>& v)                                                     { vec<T,N> out; for (int i = 0; i < N; i++) out[i] = iondrive::math::atan(v[i]); return out; }
template <class T, int N> inline vec<T,N> atan(const vec<T,N>& y,const vec<T,N>& x)                                   { vec<T,N> out; for (int i = 0; i < N; i++) out[i] = iondrive::math::atan2(y[i],x[i]); return out; }
template <class T, int N> inline vec<T,N> pow(const vec<T,N>& a,const vec<T,N>& x)                                    { vec<T,N> out; for (int i = 0; i < N; i++) out[i] = iondrive::math::pow(a[i],x[i]); return out; }
//template <class T, int N> inline vec<T,N> exp2(const vec<T,N>& v)                                                     { vec<T,N> out; for (int i = 0; i < N; i++) out[i] = exp2(v[i]); return out; }
//template <class T, int N> inline vec<T,N> log2(const vec<T,N>& v)                                                     { vec<T,N> out; for (int i = 0; i < N; i++) out[i] = log2(v[i]); return out; }
//template <class T, int N> inline vec<T,N> exp(const vec<T,N>& v)                                                      { vec<T,N> out; for (int i = 0; i < N; i++) out[i] = iondrive::math::exp(v[i]); return out; }
//template <class T, int N> inline vec<T,N> log(const vec<T,N>& v)                                                      { vec<T,N> out; for (int i = 0; i < N; i++) out[i] = iondrive::math::log(v[i]); return out; }
//template <class T, int N> inline vec<T,N> sqrt(const vec<T,N>& v)                                                     { vec<T,N> out; for (int i = 0; i < N; i++) out[i] = iondrive::math::sqrt(v[i]); return out; }
//template <class T, int N> inline vec<T,N> inversesqrt(const vec<T,N>& v)                                              { vec<T,N> out; for (int i = 0; i < N; i++) out[i] = inversesqrt(v[i]); return out; }
//template <class T, int N> inline vec<T,N> abs(const vec<T,N>& v)                                                      { vec<T,N> out; for (int i = 0; i < N; i++) out[i] = iondrive::math::fabs(v[i]); return out; }
//template <class T, int N> inline vec<T,N> sign(const vec<T,N>& v)                                                     { vec<T,N> out; for (int i = 0; i < N; i++) out[i] = sign(v[i]); return out; }
//template <class T, int N> inline vec<T,N> floor(const vec<T,N>& v)                                                    { vec<T,N> out; for (int i = 0; i < N; i++) out[i] = iondrive::math::floor(v[i]); return out; }
//template <class T, int N> inline vec<T,N> ceil(const vec<T,N>& v)                                                     { vec<T,N> out; for (int i = 0; i < N; i++) out[i] = iondrive::math::ceil(v[i]); return out; }
//template <class T, int N> inline vec<T,N> fract(const vec<T,N>& v)                                                    { vec<T,N> out; for (int i = 0; i < N; i++) out[i] = fract(v[i]); return out; }
template <class T, int N> inline vec<T,N> mod(const vec<T,N>& v, const T& t)                                          { vec<T,N> out; for (int i = 0; i < N; i++) out[i] = glsl::mod(v[i],t); return out; }
template <class T, int N> inline vec<T,N> mod(const vec<T,N>& y,const vec<T,N>& x)                                    { vec<T,N> out; for (int i = 0; i < N; i++) out[i] = glsl::mod(y[i],x[i]); return out; }
template <class T, int N> inline vec<T,N> min(const vec<T,N>& v, const T& t)                                          { vec<T,N> out; for (int i = 0; i < N; i++) out[i] = glsl::min(v[i],t); return out; }
template <class T, int N> inline vec<T,N> min(const vec<T,N>& y,const vec<T,N>& x)                                    { vec<T,N> out; for (int i = 0; i < N; i++) out[i] = glsl::min(y[i],x[i]); return out; }
template <class T, int N> inline vec<T,N> max(const vec<T,N>& v, const T& t)                                          { vec<T,N> out; for (int i = 0; i < N; i++) out[i] = glsl::max(v[i],t); return out; }
template <class T, int N> inline vec<T,N> max(const vec<T,N>& y,const vec<T,N>& x)                                    { vec<T,N> out; for (int i = 0; i < N; i++) out[i] = glsl::max(y[i],x[i]); return out; }
template <class T, int N> inline vec<T,N> clamp(const vec<T,N>& v, const T& minVal, const T& maxVal)                  { vec<T,N> out; for (int i = 0; i < N; i++) out[i] = clamp(v[i],minVal,maxVal); return out; }
template <class T, int N> inline vec<T,N> clamp(const vec<T,N>& v, const vec<T,N>& minVal, const vec<T,N>& maxVal)    { vec<T,N> out; for (int i = 0; i < N; i++) out[i] = clamp(v[i],minVal[i],maxVal[i]); return out; }
template <class T, int N> inline vec<T,N> mix(const vec<T,N>& a,const vec<T,N>& b, const T& t)                        { vec<T,N> out; for (int i = 0; i < N; i++) out[i] = mix(a[i],b[i],t); return out; }
template <class T, int N> inline vec<T,N> mix(const vec<T,N>& a,const vec<T,N>& b, const vec<T,N>& c)                 { vec<T,N> out; for (int i = 0; i < N; i++) out[i] = mix(a[i],b[i],c[i]); return out; }
template <class T, int N> inline vec<T,N> step(const vec<T,N>& edge, const vec<T,N>& x)                               { vec<T,N> out; for (int i = 0; i < N; i++) out[i] = step(edge[i],x[i]); return out; }
template <class T, int N> inline vec<T,N> smoothstep(const vec<T,N>& edge0, const vec<T,N>& edge1, const vec<T,N>& x) { vec<T,N> out; for (int i = 0; i < N; i++) out[i] = smoothstep(edge0[i],edge1[i],x[i]); return out; }

template <class T, int N> inline vec<T,N> normalize(const vec<T,N>& v)                                                { T len = length(v); if (len == 0) return v; return v/len; }
template <class T, int N> inline vec<T,N> faceforward(const vec<T,N>& n, const vec<T,N>& i, const vec<T,N>& nref)     { return dot(nref,i) < 0 ? n : -n; }
template <class T, int N> inline vec<T,N> reflect(const vec<T,N>& i, const vec<T,N>& n)                               { return i - dot(n,i)*n*T(2); }

template <class T, int N> inline T length(const vec<T,N>& v)                                                          { return T((iondrive::math::sqrt(dot(v,v)))); }
template <class T, int N> inline T distance(const vec<T,N>& a, const vec<T,N>& b)                                     { return length(a-b); }
template <class T, int N> inline T dot(const vec<T,N>& a, const vec<T,N>& b)                                          { T tmp; for (int i=0;i<N;i++) tmp += a[i]*b[i]; return tmp; }

// partial specialization for CrossProduct: it only exists in 3D vectors
template <class T> inline vec<T,3> cross(const vec<T,3>& a,const vec<T,3>& b) {
	vec<T,3> result;
	
	result[0] = a[1]*b[2] - a[2]*b[1];
	result[1] = a[2]*b[0] - a[0]*b[2];
	result[2] = a[0]*b[1] - a[1]*b[0];

	return result;
}

/* boolean vectors - there are currently no plans to try to implement
 * these boolean vectors without the use of temporaries. If speed on
 * these functions becomes an issue we'll see about that :-)
 */

// partial specialization for boolean vectors
template <class T, int N> inline vec<bool,N> lessThan(const vec<T,N>& a, const vec<T,N>& b)         { vec<bool,N> out; for (int i = 0; i < N; i++) out[i] = a[i] < b[i]; return out; }
template <class T, int N> inline vec<bool,N> lessThanEqual(const vec<T,N>& a, const vec<T,N>& b)    { vec<bool,N> out; for (int i = 0; i < N; i++) out[i] = a[i] <= b[i]; return out; }
template <class T, int N> inline vec<bool,N> greaterThan(const vec<T,N>& a, const vec<T,N>& b)      { vec<bool,N> out; for (int i = 0; i < N; i++) out[i] = a[i] > b[i]; return out; }
template <class T, int N> inline vec<bool,N> greaterThanEqual(const vec<T,N>& a, const vec<T,N>& b) { vec<bool,N> out; for (int i = 0; i < N; i++) out[i] = a[i] >= b[i]; return out; }
template <class T, int N> inline vec<bool,N> equal(const vec<T,N>& a, const vec<T,N>& b)            { vec<bool,N> out; for (int i = 0; i < N; i++) out[i] = a[i] == b[i]; return out; }
template <class T, int N> inline vec<bool,N> notEqual(const vec<T,N>& a, const vec<T,N>& b)         { vec<bool,N> out; for (int i = 0; i < N; i++) out[i] = a[i] != b[i]; return out; }

template <int N> inline bool any(const vec<bool,N>& v) { for (int i=0;i<N;i++) if ( v[i]) return true;  return false; }
template <int N> inline bool all(const vec<bool,N>& v) { for (int i=0;i<N;i++) if (!v[i]) return false; return true;  }

// C++ LIMIT: Because gcc reserves not(), so we cannot use that GLSL func. we
// use operator!() instead.
#ifndef __GNUC__
template <int N> inline vec<bool,N> not(const vec<bool,N>& v) {
	vec<bool,N> out;
	for (int i = 0; i < N; i++) out[i] = !v[i];
	return out;
}
#endif
// this operator!() can be used instead of the above not()
template <int N> inline vec<bool,N> operator!(const vec<bool,N>& v) {
	vec<bool,N> out;
	for (int i = 0; i < N; i++) out[i] = !v[i];
	return out;
}

//
// additional non-standard functions
//

template <class T,int N> inline T squared_length(const vec<T,N>& v)   { T tmp; for (int i=0;i<N;i++) tmp += v[i]*v[i]; return tmp; }
template <class T> inline vec<T,3> unit_cross(const vec<T,3>& a,const vec<T,3>& b)
{
	vec<T,3> result;

	result[0] = a[1]*b[2] - a[2]*b[1];
	result[1] = a[2]*b[0] - a[0]*b[2];
	result[2] = a[0]*b[1] - a[1]*b[0];

	result.normalize();
	
	return result;
}
// partial specialization for saru::ColorRGBA::ScaleByMax
template <class T> inline vec<T,3> scale_by_max(const vec<T,3>& v)
{
    vec<T,3> result;
    
    T max = v[0];
    if (v[1] > max) max = v[1];
    if (v[2] > max) max = v[2];
    
   for (int i=0;i<3;i++) result[i] = v[i];
 
   if (max > 1.0f)
   {
	    float invmax = 1.0f/max;
	    for (int i=0;i<3;i++) result[i] *= invmax;
   }
 
 
    return result;
}
template <class T> inline vec<T,4> scale_by_max(const vec<T,4>& v)
{
    vec<T,4> result;

    T max = v[0];
    if (v[1] > max) max = v[1];
    if (v[2] > max) max = v[2];
    
    for (int i=0;i<4;i++) result[i] = v[i];
    
    if (max > 1.0f)
    {
        float invmax = 1.0f/max;
        for (int i=0;i<3;i++) result[i] *= invmax;
    }
    
    // clamp alpha to [0,1]
    if (v[3] > 1.0f)
        result[3] = 1.0f;
    else if (v[3] < 0.0f)
        result[3] = 0.0f;
    return result;
}



/************************************************************
 *
 * matrix class
 * 
 */




template <class T, int N, int M = N>
#ifndef INCLUDED_GLSL_NO_ARRAY
class mat : public std::tr1::array<T,N*M>
#else
class mat
#endif
{
#ifndef INCLUDED_GLSL_NO_ARRAY
	//
#else
	private:		
		T data[N*M];

	public:
		// *all* access is done through this operator
		inline T& operator[](const int i)       { return data[i]; }
		inline const T&  operator[](const int i) const { return data[i]; }
#endif
		
		
	public:

		inline mat() { for (int i = 0; i < N*M; i++) (*this)[i] = T(); }

		inline T& element(const int i, const int j) { return (*this)[i*N+j]; }
		inline const T& element(const int i, const int j) const { return (*this)[i*N+j]; }
};
   
// use init() to initialize explicitly the 2,3,4D case

template <class T> inline mat<T,2> init_mat(const T& m0,
					       const T& m1,
					       const T& m2,
					       const T& m3)
{
	mat<T,2> m;

	m[0] = m0;
	m[1] = m1;
	m[2] = m2;
	m[3] = m3;

	return m;
}
template <class T> inline mat<T,3> init_mat(const T& m0,
					       const T& m1,
					       const T& m2,
					       const T& m3,
					       const T& m4,
					       const T& m5,
					       const T& m6,
					       const T& m7,
					       const T& m8)
{
	mat<T,3> m;

	m[0]=m0;
	m[1]=m1;
	m[2]=m2;
	m[3]=m3;
	m[4]=m4;
	m[5]=m5;
	m[6]=m6;
	m[7]=m7;
	m[8]=m8;

	return m;
}
template <class T> inline mat<T,4> init_mat(const T& m0,
					       const T& m1,
					       const T& m2,
					       const T& m3,
					       const T& m4,
					       const T& m5,
					       const T& m6,
					       const T& m7,
					       const T& m8,
					       const T& m9,
					       const T& m10,
					       const T& m11,
					       const T& m12,
					       const T& m13,
					       const T& m14,
					       const T& m15)
{
	mat<T,4> m;

	m[0]=m0;
	m[1]=m1;
	m[2]=m2;
	m[3]=m3;
	m[4]=m4;
	m[5]=m5;
	m[6]=m6;
	m[7]=m7;
	m[8]=m8;
	m[9]=m9;
	m[10]=m10;
	m[11]=m11;
	m[12]=m12;
	m[13]=m13;
	m[14]=m14;
	m[15]=m15;

	return m;
}

// REQUIRED: matrix operations

template <class T, int N> inline mat<T,N> operator+(mat<T,N> a, mat<T,N> b) { mat<T,N> out; for (int i = 0; i < N*N; i++) out[i] = a[i]+b[i]; return out; }
template <class T, int N> inline mat<T,N> operator-(mat<T,N> a, mat<T,N> b) { mat<T,N> out; for (int i = 0; i < N*N; i++) out[i] = a[i]-b[i]; return out; }
template <class T, int N> inline mat<T,N> operator+(const T& a, mat<T,N> m) { mat<T,N> out; for (int i = 0; i < N*N; i++) out[i] = a+m[i]; return out; }
template <class T, int N> inline mat<T,N> operator+(mat<T,N> m, const T& a) { return a+m; }
template <class T, int N> inline mat<T,N> operator-(const T& a, mat<T,N> m) { mat<T,N> out; for (int i = 0; i < N*N; i++) out[i] = a-m[i]; return out; }
template <class T, int N> inline mat<T,N> operator-(mat<T,N> m, const T& a) { mat<T,N> out; for (int i = 0; i < N*N; i++) out[i] = m[i]-a; return out; }
template <class T, int N> inline mat<T,N> operator*(const T& a, mat<T,N> m) { mat<T,N> out; for (int i = 0; i < N*N; i++) out[i] = a*m[i]; return out; }
template <class T, int N> inline mat<T,N> operator*(mat<T,N> m, const T& a) { return a*m; }

// matrix-vector multiplication
template <class T, int N> inline vec<T,N> operator*(mat<T,N> m, const vec<T,N>& v)
{
	vec<T,N> out;

	for (int j = 0; j < N; j++)
	{
		out[j] = T();
		
		for (int i = 0; i < N; i++)
			out[j] += m[j+N*i] * v[i];
	}
	
	return out;
}

// matrix-matrix multiplication
template <class T, int N> inline mat<T,N> operator*(mat<T,N> a, mat<T,N> b)
{
	mat<T,N> out;

	for (int k = 0; k < N; k++)
	{
		for (int j = 0; j < N; j++)
		{
			out[k*N+j] = 0;
			
			for (int i = 0; i < N; i++)
				out[k*N+j] += a[j*N+i] * b[k*N+i];
		}
	}
	
	return out;
}

// additional non-square matrices multiplications
/*
 * Here are two multiplications I really *need* to see implemented in
 * the standard header. Most simulations need a list of array values
 * for time simulations, for example.
 */

// vec_100 = mat_100_200 * vec_200 type multiplication
template <class T, int N, int M>
inline glsl::vec<T,N>
operator*(glsl::mat<T,N,M> m, glsl::vec<T,M> v)
{
	return glsl::vec<T,N>();
}
// mat_10_20 = mat_10_15 * mat_15_20 type multiplication 
template <class T, int N, int M, int O>
inline glsl::mat<T,N,O>
operator*(glsl::mat<T,N,M> A, glsl::mat<T,M,O> B)
{
	return glsl::mat<T,N,O>();
}




// builtins (there is only one of this in original GLSL)
template <class T, int N, int M> inline mat<T,N,M> matrixcompmult(mat<T,N,M> a, mat<T,N,M> b) { mat<T,N,M> out; for (int i = 0; i < N*M; i++) out[i] = a[i]*b[i]; return out; }

//
// additional non-standard functions
//

template <class T, int N, int M> inline mat<T,M,N> transpose(mat<T,N,M> m)
{
	mat<T,M,N> out;

	for (int j = 0; j < M; j++)
		for (int i = 0; i < N; i++)
			out[j*M+i] = m[i*N+j];
	
	return out;
}

template <class T, int N> inline mat<T,N> diag(const vec<T,N>& v)
{
	mat<T,N> out;			// default ctor sets all to 0
	
	for (int i = 0; i < N; i++)
		out[i*N+i] = v[i];
	
	return out;
}
template <class T, int N> inline mat<T,N> diag(const T& v, const int n=N)
{
	mat<T,N> out;			// mat3 foo = diag(1,3);=makeidenity; foo=diag(0,3); == makezero

	for (int i = 0; i < N; i++)
		out[i*N+i] = v[i];

	return out;
}




/************************************************************
 *
 * These are all the standard GLSL types as typedefs
 * 
 */




typedef vec<float,2>  vec2;
typedef vec<int,  2> ivec2;
typedef vec<bool, 2> bvec2;

typedef vec<float,3>  vec3;
typedef vec<int,  3> ivec3;
typedef vec<bool, 3> bvec3;

typedef vec<float,4>  vec4;
typedef vec<int,  4> ivec4;
typedef vec<bool, 4> bvec4;

typedef mat<float,2>  mat2;
typedef mat<float,3>  mat3;
typedef mat<float,4>  mat4;


/************************************************************
 *
 * The following part just handles text output for streams
 *
 */

template<class T, int N> std::ostream& operator<<(std::ostream& os, vec<T,N> v)
{
	os << '(' << v[0];
	for (int i=1;i<N;i++) os << ',' << v[i];
	os << ')';

	return os;
}
template<class T, int N, int M> std::ostream& operator<<(std::ostream& os, mat<T,N,M> A)
{
	os << '(';
	
	os << '(' << A.element(0,0);
	for (int i=1;i<M;i++) os << ',' << A.element(0,i);
	os << ')';
	
	for (int rows = 1; rows < N; rows++)
	{
		os << ',' << ' ';
		
		os << '(' << A.element(rows,0);
		for (int i=1;i<M;i++) os << ',' << A.element(rows,i);
		os << ')';
	}
	

	os << ')';
	
	return os;
}



} // end namespace glsl
} // end namespace iondrive

#endif /* INCLUDED_GLSL_H */



